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ABSTRACT 

We propose a robust ensemble filtering scheme based on the filtering theory. The optimal 
ifoo filter is derived by minimizing the supremum (or maximum) of a predefined cost function, 
a criterion different from the minimum variance used in the Kalman filter. By design, the 
filter is more robust than the Kalman filter, in the sense that the estimation error in the 
ifoo filter in general has a finite growth rate with respect to the uncertainties in assimilation, 
except for a special case that corresponds to the Kalman filter. 

The original form of the filter contains global constraints in time, which may be 
inconvenient for sequential data assimilation problems. Therefore we introduce a variant 
that solves some time-local constraints instead, and hence we call it the time-local filter 
(TLHF). By analogy to the ensemble Kalman filter (EnKF), we also propose the concept 
of ensemble time-local ifoo filter (EnTLHF). We outline the general form of the EnTLHF, 
and discuss some of its special cases. In particular, we show that an EnKF with certain 
covariance inflation is essentially an EnTLHF. In this sense, the EnTLHF provides a general 
framework for conducting covariance inflation in the EnKF-based methods. We use some 
numerical examples to assess the relative robustness of the TLHF/EnTLHF in comparison 
with the corresponding KF/EnKF method. 



1. Introduction 



The Kalman filter (KF) (jKalmanlll960l ) is a sequential data assimilation algorithm. For 



linear stochastic systems, it can be shown that the K 



minimizes the variance of the estimation error fjSimon 



is an optimal linear estimator that 
20061 . ch. 5). Because of its relative 



simplicity in implementation, the KF is suitable for many data assimilation problems. How- 
ever, for high dimensional systems such as weather forecasting models, direct application 
of the KF is prohibitively expensive as it involves manipulating covariance matrices of the 
system states. For this reason, different modifications of the KF were proposed to reduce 



the computational cost. These include various ensemble Kalman filters (EnKF) (Anderson 



2001[ iBishop et al.ll2001t iBureers et al.lll998uEvensedll994|; | Evensen and van Leeuwenlll996 ; 



Houtekamer and Mitchell 



1998 



2007 



ters ( Cohn and Tod ling 



1998 



1996 



Verlaan and Heemink 



Zupanski 



Whitaker and Hamill 



H 



oteit et al. 



2001 



2002 



20021). the error subsp a ce-based fil 



Luo and Moroz 



19971 ) , and filters based on other strategies (IBeezley and Mandel 



2009 



Pham et al 



20051 ) . to name but a few. A detailed 



yond the scop e of th is work. Readers are referred to 



Tippett et al 



descript i on of t he above filter s is be 



Evensenl (120031 ) 



Nerger et al. 



(120051) 



(120031 ) for reviews of some of the aforementioned filters. Roughly speaking, 



these modifications exploit the information of a subset in the state space of a dynamical 
system, while the information of the complement set is considered less influential, and thus 
ignored. Consequently, the computations of these modified filters are normally conducted 
on the chosen subsets, instead of the whole state space, so that their computational costs 
are reduced. For simplicity, we may sometimes abuse the terminology by referring to all the 
aforementioned filters as the EnKF-based methods (EnKF methods for short). 

The KF and the EnKF are among the family of Bayesian filters that adopt Bayes' rule 
to update background statistics to their analysis counterparts. In these filters, one needs 
to make certain assumptions on the statistical properties (e.g., pdfs or moments) of both 
the dynamical and observation systems. In reality, however, these assumptions may not be 



accurate, so that a Bayesian filter may fail to achieve good performance with mis-specified 



statistical information f jSchlee et al.l 119671 ). For example, if implemented straightforwardly, 
an EnKF with a relatively small ensemble s ize may produce inaccurate estimations of covari- 



ance matrices (IWhitaker and Hamill 



20021 ). This could degrade filter performance, or even 



cause filter divergence. As a remedy, in practi ce it is customary to conduct c ovariance in 



flat ion anc. 



2001 



2009 



ocalization to relie ve this problem ( lAnderson and Anderson 



1999 



Hamill et al 



Van Leeuwen 



20091 ). 



In contrast, robust filters emphasize on the robustness of their estimates, so that they may 



he estimation strategies o f 



1987 



Simon 



2006|) 



have better tolerances to possible uncertainties in assimilation, 
robust filters are different from Bayes' rule. Take the filter (jFrancis 
one of the robust filters, as an example. The filter (HF) does not assume to exactly know 
the statistical properties of a system in assimilation. Instead, it accepts the possibility that 
one may only have incomplete information of the system. Consequently, rather than looking 
for the best possible estimates base d on Bayes' ru le, the optimal filter employs a robust 



strategy, namely, the minimax rule ( jBurger 



1985 



ch. 4), to update its background statistics. 



This robustness may be of interest in practical situations. For example, for data assimilation 
in earth systems, the system models are often not the exact descriptions of the underlying 
physical proce sses, and it is challen ging to characterize the properties of the corresponding 

and the references therein). Given an imperfect model, 



model errors (jWang and Cai 



2008 



the estimation error of the HF in general grows with the uncertainties in assimilation at a 
finitely bounded rate (except for the special case when the HF reduces to the KF itself), 
while the estimation error of the KF does not possess such a guarantee. 

In this work we propose a variant of the HF, called the time-local HF (TLHF), to avoid 
solving global constraints as in the HF. By analogy to the EnKF, we further introduce the 



ensemble TLHF (EnTLHF) for data assimilation in high dimensional systems. We show 
that the EnTLHF can be constructed based on the EnKF, and thus the computational 
complexity of the EnTLHF is in general comparable to that of the EnKF. We also show 
that some specific forms of the EnTLHF have connections with some EnKFs equipped with 
certain covariance inflation techniques. More generally, we show that an EnKF with a certain 
covariance inflation technique is in fact an EnTLHF. 

The organization of this work is as follows. Section [2] presents the data assimilation 
problem and discusses its solutions in terms of the KF and the HF, respectively. Section [3] 
introduces the TLHF as a variant of the HF, and its ensemble form, the EnTLHF. Section 
H] discusses some specific forms of the EnTLHF and shows their connections with some of 
the EnKF methods with covariance inflation. In section [51 we use some numerical examples 
to show the relative robustness of the TLHF/EnTLHF in comparison to the corresponding 
KF/EnKF method. 

2. Problem statement 

We consider the state estimation problem in the following scenario: 



Xj = M it i-i (Xi_i) + Ui , (la) 

y, = T-ii (xj) + Vi , (lb) 

E Ui = 0; E( Ui uJ) = , (lc) 

E Vi = 0; E( Vi vJ) =^Ri, (Id) 

E(u,vj)=0 Vz,j. (le) 
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Eqs. ([Taj) and (llb|) represent the m x - dimensional dynamical system and the m y - dimensional 
observation system, respectively, where x, denotes the m^-dimensional state vector, the 
corresponding m^-dimensional observation, Aii t i+i : IR" 1 * — >• M. mx the transition operator, 
and Hi : M. mx — >• ~R my the observation operator, at time %. For convenience, in this section we 
assume the systems in Eqs. ( fTaj) and (llbl) are linear so that A^i.i+i and Hi are two matrices. 
The estimation problem in the presence of nonlinearity will be addressed in terms of the 
EnTLHF, as will be discussed later. 

Eqs. (ITc|) - (fTe|) imply that the m^-dimensional dynamical noise Uj and the m^-dimensional 
observation noise are uncorrelated white noisqj, with zero mean and covariances Qj and 
R;, respectively. 5y denotes the Kronecker delta function such that 5ij = 1 if i = j, Sy — 
otherwise. 

In what follows, we discuss two filtering approaches as the solutions of the above state 
estimation problem: The KF which is based on the criterion of minimizing the variance of 
the estim ation error (equi valent to applying Bayes' rule to update background statistics as 



shown in 



Jazwinski 



1970 



ch. 7), and the optimal HF which is based on the criterion of 

minimizing the supremum (or maximum) of the ratio of the "energy" of the estimation error 

to the "energy" of the uncertainties in data assimilation (to be made clear shortly). In what 

follows we outline the main results of the KF and the HF. For more details, readers will be 

referred to appropriate references. 

1 The deducti on will be sim ilar in case that and are correlated colored noise. Readers are referred 
to, for example, ISimonl (|2006l ch. 7) for the details. 
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a. Kalman filter as a solution 



The KF algo rithm involves prediction and filtering steps, whose deductions can be found 



in, for example, 



Simon! ( 12006 



ch. 3). When the KF is applied to assimilate the system in 



Eq. (JTJ), these steps follows. 



Prediction step: Propagate the analysis x°_ 1 at the (i — l)th cycle and the associated 
analysis error covariance Pf_ x forward to obtain the background at the next cycle 
(Eq. ( |2al) ) and the associated background error covariance P^ 1 (Eq. ( |2bl) ). respectively. 



P^A^wPSLi^w-i + Qi 



(2a) 
(2b) 



Filtering step: With a new incoming observation y^, update the background x^ and 
the associated error covariance P^ to their analysis counterparts, x" (Eq. ( 15a)) ) and P® 
(Eq. (155)) ). respectively. 



x? = x* + K, (y, - ffcx}) , 
P° = P* - K^P^ , 

K, = p^n^p^f +p,) _i 



(3a) 
(3b) 
(3c) 



where K 7 in Eq. (15cl) is the Kalman gain. Alternatively, Eqs. ( I3bj) and ( 15c)) can be reformu- 



lated as fjSimon 



20061 . ch. 6) 



>a\ — 1 



,6\-l 



("Hj) T (Ri) 1 Hj 



K, = P^f(R,)- 1 



(4a) 
(4b) 



Eq. (j3J) will be used to simplify the presentation of the HF later on. 

To better understand the difference between the KF and the HF, we extend our discussion 
to a slightly more general situation. Instead of looking for an estimate of the system state 
Xj, we are interested in estimating some m z -dimensional linear transform z, of Xj, in terms 
of 

Zj = LjXj , (5) 

with Lj being a user-defined m z x m x matrix. In Eq. (jSJ), if Lj is the m x x m x identity 
matrix, then Zj reduces to Xj itself. But Lj can be something else in general. For instance, 
one may let Lj = [1, 0, ■ ■ • , 0] such that Zj is equivalent to the first element of Xj. 

In the KF, one can solve the estimation problem by minimizing the following cost function 

N N 

J? F (z%, ...,z%) = J2 j™ = E E ii z * - z "\\i ( 6 ) 

i=0 i=0 

with respect to the variables Zq, • • • , z^. Here, z" represents the analysis of the truth z,, 
while Jj^ F = E||zj — z°|||, the cost function local to the ith cycle with respect to Zj, is the 
variance of the estimation error (z, — z") at the ith assimilation cycle, with E being the 
expectation operator, and || • ||a denoting the L 2 -norm. Given an m 2 -dimensional vector 
Z = (z ± , ■ ■ ■ ,z mz ) T , ||z|| 2 = Vz r z = JJ2 Z ]- 

The estimates of the KF are obtained by sequentially minimizing the local cost function 
J^f and then propagating the resulting analysis forward to yield the background at the next 
cycle. By treating z$ = LjXj as a system state, direct sequential minimization of the local 
cost functions Jff is equivalent to running the KF recursively to assimilate the following 
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modified systems 

Zj = (L^i.i-iLr 1 ) Zj_! + LiUj , (7a) 
y, = (T^Lr 1 ) Zi + Vi , (7b) 

where Zj and LjUj are the system state and dynamical noise at instant i, respectively, as- 
sociated with the modified transition operator LjA^j-iL" 1 , where L" 1 may have to be 
treated as the pseudo-inverse of Lj in some circumstances, y$ and Vj still correspond to the 
observation vector and noise, respectively, but the observation operator becomes T-LjhJ . 

Alternatively, one may solve the estimation problem in an indirect way, but without 
changing the systems in assimilation. To this end, one specifies a new cost function 

N N 

j« f w, ■ ■ ■ , x*) = E J S = E E H x * - x "H* ( § ) 

i=0 j=0 

with respect to the variables xg, ■ ■ • , x^r to determine the optimal estimations by sequentially 
minimizing the local cost functions J^f (equivalent to running the KF recursively to assim- 
ilate the system in Eq. (JT])), and then applies Eq. (JSJ) to obtain the estimations zg, • • ■ ,z^ 
of interest. 

The user-defined matrix Lj does not appear in Eq. (JSJ), which implies that in indirect 
estimation, the estimates x" of the KF do not depend on the choice of Lj, so that by putting 

z ■ = L,x? , (9) 

direct and indirect methods yield the same estimates z". However, for the HF to be presented 
below, we will see that the choice of Lj does affect the estimates x" in indirect estimation, 
which is clearly different from the KF. 
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b. filter as a solution 



The HF (lSimonll2006l . ch. 11 and the references therein) aims to provide robust, but not 
necessarily best, estimates. The main idea is to first recognize that in the system Eq. (JT]), 
there are three possible sources that contribute to the uncertainties in data assimilation, 
namely, the uncertainties in the initial conditions, the model error and the observation 
error. Accordingly, during an assimilation time window [0, N], these uncertainty sources are 

N 

characterized by three uncertainty "energy" terms, defined as ||x — x |R _ 1; ||uj||^_ 1 , and 

N 

S II v «IIr-15 respectively. For consistency, here we assume that an assimilation cycle starts 
from propagating an analysis forward through a dynamical system, so that the dynamical 
noise uq is included in our discussion. Given a symmetric, positive semi-definite matrix A, 



II • II a represents the weighted L 2 -norm so that ||z||a = V z T Az for a column vector z. 

Since the minimum variance criterion in the KF is consistent with Bayes's rule, it is cus- 
tomary to interpret the matrices Ao, Qi, and R« as the covariances (or their estimates) of 
the corresponding random vectors in probabilistic language. In the context of the HF, how- 
ever, the minimax rule is adopt ed instead, w hich does not have probabilistic interpretations 



for A , Qi, and R» in general (jSimon 



20061 ). Therefore, for conceptual distinction, in the 
context of the HF we call A , Qi, and Rj the uncertainty weight matrices (uncertainty ma- 
trices for short) with respect to the initial conditions, the model error, and the observation 
error, respectively, and their inverses, Aq , Q^ 1 , and R~\ the corresponding information 
matrices. The uncertainty matrices A , Qi, and R, can be user-defined, although in practice 
it is convenient to let them be equal to their covariance counterparts (or the estimates) in the 
KF. Under this choice, the difference between the estimates of the KF and the HF becomes 
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clear, as will be shown later. 

To provide robust estimates, the HF requires that the total "energy" of the estimation 



N 



errors, in terms of ^ ||zj — zf |||., be no larger than the sum of the above three uncertainty 



i=0 



"energy" terms times a constant factor I/7, i.e., 



N 



1 



/V 



A 



hi - tfWl < i- [ llxo - xolli-i + 2J IKIIq- + U llvill^ I • 

i=0 ^ \ i=0 i=0 

In (fTOl) . the weight ma trix Sj is als o user-chosen, which can be used to guide the filter's 



behavior in estimation (I Simon 



20061 . ch. li; 



To solve ( II 01) . we first define the following cost function J. 



HF 



A 



J 



HF 



XI \\ Z i Z i\\si 
i=0 



AT N 

Ixo-xolli^ + EKII^ + I^NI^ 

j—Q I j 



(Hi 



which is a function of the variables x , {uj}, {v^} and {z"} q Then the inequality (fTUj) is 
equivalent to J^ F < I/7, except when the uncertainty "energy" term 



A? 



A' 



|X " Xo||^-i + ii^Hq- 1 + 1] H V 



i=0 i=0 

is zero. In this exceptional case, there still exists a HF solution to make the inequality (TTOj) 
hold. For example, if A , Q« and Rj are all +00 so that their inverses are zero, we can 
choose S» = so that the left hand side (lhs) of (fTUj) is also zero (for consistency, we define 
the products oo-0 = 0- oo = 0to avoid possible problems). Similarly, if xo — xo, u« and Vj 
are all zero, we have the ideal solution z" = Zj to make the lhs of (llOp become zero. 
Let 7* be the value such that 



— = inf sup Jf F , i< N, 

1* {<} x ,{ Ul },{v s } 



'12) 



2 Xo is the prior knowledge of xo, which is assumed to be available, and therefore not influenced by the 
filter designer. 



10 



where sup J^ F calculates the supremum of the cost function Jf F with respect to the 

X ,{Uj},{v;} 

variables x , {u^}, {vj} (which is a function of {z"}), and inf evaluates the infimum with 
respect to {z"} along the previously obtained supremum plane of Xo, {u^}, {v.;}. We say that 
the optimal HF is achieved if 7 = 7*. In this sense, the evaluation of 7* is an application of 
the minimax rule, a st rategy that a ims to provide robust estimates and is different from its 



Bayesian counterpart (IBurger 



1985 



% chs. 4 and 5). 

In practice, it is difficult to evaluate the exact value of 7*, since, by the definition in 
Eq. (TT2|) . 7* depends not only on the initial conditions and the dynamical and observation 
systems, but also on the length iV of the assimilation time window. A more convenient 
strategy is to choose a value 7 satisfying I/7* < I/7 < +00, so that it gu arantees that 



2006 



there exists a (sub-optimal) HF solution {z"} satisfying the inequality (fTOl) (ISimon 
ch. 11). In contrast, if I/7 < I/7*, then there is no guarantee (although it is not impossible) 
that there exists such a HF solution. For instance, in the extreme event with the variables 
xo, {uj}, {vj} located on their supremum plane, then by the definition in Eq. (fl2l) . there is 
no HF solution {z"} to make J^ F < I/7 < I/7*. For reference, hereafter we call 7 the 
performance level of the HF. 

The inequality ffTOl) can be solved through dynamic constrained optimization, with Eqs. (fTa| 



and ( llbl) being t he co nstraints at different assimilation cycles. For details, readers are re- 



ferred to 



Simonl (120061 . ch. 11). For convenience of comparison, we also split the algorithm 



into prediction and filtering steps. 



Prediction step: As in the KF, we also propagate the analysis forward to produce the 
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background at the next cycle. 

x£ = M M -ixti, (13a) 
A^Mi^AUMj^ + Qi. (13b) 

Filtering step: With a new observation y i; we update the background to the analysis. 

^ = ^ + G l (y l -'H l ^) , (14a) 

(A?)- 1 = (A, 6 )" 1 + (^m)- 1 ^ - iLfS^ , (14b) 

G t = A^HR,)- 1 , (14c) 
subject to the constraints 

(A?)" 1 = (A^)" 1 + (n i ) T (R i )- 1 n t - TLfSiL* > . (15) 

Here A, denotes the uncertainty matrix, analogous to the covariance matrix Pj in the KF, 
and Gj represents the gain matrix, analogous to the Kalman gain Kj. The symbol "> 0" 
means that (A")" 1 must be positive semi-definite. After obtaining x", one can apply Eq. (j9]) 
to estimate Zj. 

Comparing Eq. ffTB"]) with Eq. (J2]) , one can see that the prediction steps of the KF and the 
HF are the same. Furthermore, the update formula Eq. (j!4ap of the HF is a linear estimator 
as in the KF, but in general with a different gain matrix Gj. However, if 7 = in Eq. ( 114b|) . 



then Eqs. (I14a[) . ( 114b|) and f ll4cj) reduce to Eqs. (l3aj) . (|4aj) and ( libj) . respectively. In this 



case, the HF is equivalent to the KF. Therefore, the KF can be considered as a special case 
of the HF with a performance level 7 = 0. 

A further examination of Eq. ( 114]) reveals that the user-defined matrix Lj is involved in 
the evaluation of A", hence those of Gj and x". This implies that, unlike the KF, the choice 
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of Lj affects the estimates of the HF. In addition, let 

(S«)- 1 = (A l 6 )- 1 + (K i f(R i )- 1 ^, (16) 

then following Eq. ( Hal) , it is clear that is an uncertainty matrix obtained by updating 
A^ through the KF, while the information matrix 

(An^ 1 = (Sr)^-7LfS l L l <(S«)- 1 . (17) 

Therefore, the HF appears more conservative, in the sense that it tends to make its analysis 
uncertainty matrix A" larger than its counterpart in the KF, given the same background 
uncertainty matrix A^. Consequently, when using Eq. (I14a|) to update the background 
to the analysis, the HF allocates more weight to the observation y, than the KF does, 
which may be preferred when the background is not very reliable because of the uncounted 



sources of uncertainties (similar arguments can also be found in, for example, 



1970 : 



Van Leeuwen 



Jazwinski 



20091 ) □. In fact, as will be shown later, this conservativeness exhibits 
connections to covariance inflation techniques adopted in some EnKF methods. 

The presence of the term — ^yhj SjLj in Eq. (114bl) reflects the fundamental difference 
in the estimation strategies employed in the KF and the HF. As pointed out previously, 
the KF assumes that one has sufficient information of the statistical properties of both 
the dynamical and observation systems. Under this assumption, the KF looks for the best 
possible estimates under a certain optimality criterion (e.g., minimum variance or maximum 
likelihood). Therefore, the KF can achieve good performance if the statistical description of 



3 If, in contrast, the observation is very unreliable, then one may choose a negative value for 7 such that 
the background has relatively more weight in the update. In this work we confine ourselves to the scenario 
7 > 0. 

13 



the assimilated system is sufficiently close to the reality. However, if there exists substantial 



deviation, then the KF may have poor performance, or even diverge (see ISchlee et al.lll967 



for such examples). In contrast, the HF is more conservative than the KF since it only 
aims to provide robust, rather than best, estimates. The HF may perform worse than the 
KF if the statistical properties (e.g., pdf or moments) of the assimilated system are well 
characterized and close to the tr uth. Howeve r, if there exist more uncertainties, the HF can 



perform better than the KF (cf. 



Simon 



20061 . p. 358 for some examples). 



A further issue that may be of interest in practice is the choice of the term — 7L? 1 SjLj, 
such that one may improve the performance of the filter under other measures, e.g., root 
mean squared errors (RMSE). This is possible, since the HF only accounts for the robustness 
property, and one has certain freedom in choosing the forms of the quantities 7 and Sj □, 
as long as — 7L^SjLj satisfies the constraints in the HF. In other words, the robustness 
requirement in general does not yield unique choices of 7 and Sj. One may impose other 
objective functions to the HF to improve its performance un der the relev ant measures. One 



such example is the mixed (nonlinear) Kalman and Hoc filter ( iSimon 



2006 



ch. 12), where one 



tries to minimize a least square cost function as in the KF, while imposing certain robustness 
on the estimates in th e sense of Eq. f fTOj) . Other examples can be found in Anderson's works 



([Anderson 



2007 



20091 ) on adaptive covariance inflation in the EnKF, where the criterion that 



the pdf of the covariance inflation factor be maximized is used. This point will become clear 
after we establish the connection between robustness and covariance inflation in subsequent 
sections. 



'Here we assume Lj is determined by the practical need given in Eq. ([5]). 
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3. Time-local ensemble filter 

The HF has to satisfy the inequality constraints in ( fl5l) . which makes it challenging and 
inefficient for sequential data assimilation in certain circumstances. To see this, suppose that 
for % — 0, • • • ,N, the HF has an admissible solution {x"}, with a pre-specified performance 
level 7 satisfying all the inequality constraints in ( fi"5l) . When extending the time horizon 
from N to N + 1, 7 may not satisfy the constraint at N + 1. As a result, one has to choose a 
smaller value for 7 and re-start the assimilation in the new time window [0, N + 1], resulting 
in a different filter solution for i — 0, • • • , N. 

Alternatively, one may keep the solution between % = 0, • • • ,N unchanged. From N + 1, 
one uses a smaller value 7' for estimation as long as it satisfies (115j) . Once 7' violates the 
constraint for a larger iV, one adopts an even smaller performance level 7", but still keeps 
the previously obtained estimates, and so on. In what follows, we extend this idea further. 
We propose a variant of the HF, called the time- local Hqo filter (TLHF), in which we impose 
a local cost function and adopt a local performance level 7$ to solve a local constraint at 
each assimilation cycle. 

a. Time-local filter for linear systems 

We first introduce the TLHF for linear systems. The extension to nonlinear systems, 
analogous to the EnKF methods, will be presented in the next section. 
In the TLHF, we define a local cost function 




x + IKII 2 ^ + IN|j_ 



-1 



(18) 



15 



with respect to the variables Xj, u i; v, and z". Here is the background and is fixed 
relative to the ith cycle (xq = x ). For each i, J^f can be treated as a special case of the 
cost function J^ F in Eq. (TTTj) . as if N = i = 0. 
Analogous to (ITU]) , we require 

,o||2 ^ 1 ( ||„ „6||2 



K ^Xi-xJH^ + 11^11^ + 11^11^ , (19) 
where 7, is a suitable local performance level, which satisfies 

i>l^inf sup Jff, (20) 

h U Z i Xi,Uj,Vj 

with I/7* being the minimax point of J^f . In addition, 7$ also has to satisfy a local 
inequality constraint at time instant i (to be shown later). Then for all such 7^ (0 < i < N), 



we have 

N N 

k - < E - x * - ilk')- 1 + ii^Hqt 1 + INI ^ 

i=0 i=0 ,l \ 

^ / JV V AT \ 

^ m f x {-> E n - ^nk- 1 + E ii^Hq- + E h^Hr- • 

h \i=0 i=0 i=0 / 

This shows that the growth rate of the total "energy" of the estimation error is finite unless 
7i = for some i. Thus, the corresponding estimates {x"} may also provide a robust solution 
in the same sense as in the HF. 

Eq. (I2TI) bears a similar form to ([TO]) , but also exhibits a clear difference. That is, in 
Eq. (I2TI) . the total "energy" of the uncertainties includes the contribution from the un- 
certainty in specifying the background at each assimilation cycle. In contrast, in ( fTPl the 
counterpart term only represents the contribution from the uncertainty in specifying the 
initial conditions. The extra terms in Eq. (12T|) provide a possibility to take into account the 
effect (s) of nonlinearity and/or other mechanisms that contribute to the estimation errors 
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in the background, so that one does not have to significantly change the structure of the HF 
when extending it from linear systems to nonlinear ones. For example, in the presence of 
nonlinearity, there may exist extra uncertainties incurred by the chosen data assimilation al- 
gorithm itself (called algorithm uncertainty hereafter), including the linearization error when 
one uses the extended Kalman filter (EKF) to assimilate a nonlinear system, and more gen- 
erally, the approximation error when one adopts an approximation scheme in assimilation, 
such as the sampling error in the EnKF, or the rank deficiency in a reduced rank filter. These 
potential uncertainties influence the estimations of the system states, but conceptually they 
might not belong to the uncertainties in specifying the dynamical or observation systems. 
Instead, one may treat them as the uncertainties in specifying the background, an extension 
of the uncertainties in specifying the initial conditions. With this treatment, one may apply 
the TLHF to a nonlinear system in the same way as it is applied to a linear system, while 
including the uncertainties due to the effect (s) of nonlinearity and/or any other error sources 
into the category of uncertainty in sp ecifying the b ackground^ . 



Following the same deductions in 
TLHF as follows. 
Prediction step: 



Simon 



( 2006 



ch. 11), one can derive the steps of the 



(22a) 
(22b) 



5 Like the extended KF, there also exists the exte nded HF containing more thorough treatment of non- 
linearity (see, for example. IShaked and Bermanlll995l ). whose implementation, however, involves the deriva- 



tive^) of nonlinear functions and more sophisticated inequality constraints. 
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Filtering step: 

x? = x* + Gi (y, - Urf) , (23a) 

(A^ 1 = (A-) _1 + (^) T (Ri)-'^ - 7iLf S.L, , (23b) 

G, = A^CH^)- 1 , (23c) 

subject to the constraint 

(A?)" 1 = (A^)- 1 + (H t f (R,)- 1 ^ - jiLjSiLi > . (24) 

Thus compared with the HF, the TLHF only replaces the (global) performance level 7 
with the local one 7, (i — 1, • • • , A r ), without changing anything else. 

b. Ensemble time-local filter 

The ensemble time-local Hqo filter (EnTLHF) is a straightforward analogue to the EnKF 
methods. Here the principal idea is that, at the prediction step, one uses the background 
ensemble, which is the propagation of the analysis ensemble from the previous cycle, to 
estimate the background and the associated uncertainty matrix. Then, one updates the 
background uncertainty matrix to the analysis one based on an EnKF method, calculates 
the corresponding Gain matrix of the EnTLHF, and computes the analysis mean and the 
associated uncertainty matrix (cf. Eq. (128]) below). 

Concretely, let X* = {x% • : x| • = A^j i j_i(x"_ 1 •), j = 1, ••• , n} be the n-member 
background ensemble at time instant i, which is the prediction of the analysis ensemble 
X"_! = {x-i-ij,j = 1, • " " , n } & t the previous cycle, with the transition operator 
possibly being nonlinear. 
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At the prediction step of the EnTLHF: 

it b { = mean(X^) , (25a) 
A^ = cov(X^) + Q i , (25b) 

where and A.^ denote the estimations of the background and the associated uncertainty 
matrix, respectively, and may be obtained through different ways in various EnKF methods. 

At the filtering step, one uses A^ in place of the background sample covariance in an 
EnKF method, and applies the corresponding formula in the EnKF method to update A^ 
to an uncertainty matrix E". Then by Eqs. ( fl6l) and ( fTTl) . one has 

A? = (U - (Lf S.L,))" 1 £° , (26) 

where I mx denotes the m x - dimensional identity matrix. Thus in general A" can be evaluated 
based on the analysis update through an EnKF method. 

Substituting Eq. (1261) into Eq. (I23cl) . we have the estimated gain matrix 

G, = A^HR;)- 1 = (\ mx - 7l S?(Lf S.L,))" 1 K, , (27) 

where K; = ^Uj (R) -1 is the Kalman gain evaluated by an EnKF method, due to Eq. fjibl . 
Note that for simplicity, in the above deduction we have assumed that the observation 
operator Hi is linear. In the case that Hi is nonlinear, one could also derive an approximate 
formula of Eq. (127]) (e.g., by linearizing or using a Monte Carlo approach). 

In particular, if one chooses 7« = for i = 0, • • ■ , N in Eqs. ( 1261) and ( 1271) . then it is 
clear that A" and Gj are reduced to and K,;, their counterparts in the EnKF method, 
respectively. In this case, the EnTLHF is equivalent to the EnKF method without covariance 
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inflation. If ji > 0, the EnTLHF has connections to some EnKF methods with certain 
covariance inflation techniques, as will be discussed in §U 



For summary, at the filtering step of the EnTLHF: 



[S?,Ki] =EnKF(X J 6 ,Q i) ?{ i 



Gi = [l mx - 7 i S?(L i SjL 



K, 



x? = x$ + & (y 4 - 
A" = (i^ - 7iS?(Lf S,L 4 ) ) 



(28a) 
(28b) 
(28c) 
(28d) 



subject to the constraint 



A J a )- 1 = (S«)" 1 - 7l L/S l L l >0 



(29) 



where Eq. (128aj) means that X" and K, are obtained through an EnKF method, given the 
information X^ 1 , Q», and 



After obtaining x" and A", one can draw an analysis e nsemble o 



preserves the estimate x- and the uncertainty m atrix A" (jAnderson 



the system states, which 



2001 



Hoteit et al. 



2002 



Whitaker and Hamill 



Bishop et al. 



2001 



20021 ). This analysis ensemble is then integrated 



forward to start a new assimilation cycle. 



4. Some specific forms and their connections to covari- 



ance inflation 



Here we show that some specific forms of the EnTLHF exhibit connections to some 
existing EnKF methods with covariance inflation. We again assume that the observation 
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operator 7ii is linear. This choice is made only for convenience in presenting the term 
('Hj) T (Ri)~ 1 'Hi below. The same results can also be obtained if is nonlinear. 
From Eq. ( |29l) . we have 

(A?)" 1 = (A*)" 1 + (^) T (R i )- 1 ^ - 7i Lf SiLi = (S^)- 1 - 7jLf SjLj > , (30) 

where the local performance level 7$ and the information matrix Sj can be chosen by the 
designer. We consider some specific choices of the term 7«Lf SjLj, and derive the correspond- 
ing relations between Af and Sf (or A£). By treating A? as the inflated covariance of E", 
the EnTLHF can be interpreted as an EnKF method equipped with a specific uncertainty 
inflation technique. 

a. Case ^Lf S^L, = c(A J 6 )" 1 for < c < 1 
This implies that 

S, = — (LiLf )" 1 L i (A-)~ 1 (Lf )(LjLf ) _1 , (31) 

H 

and 

(A^)" 1 = (1 - c)(K\Y l + (K i ) T (R i )- 1 ^ . (32) 

In particular, if Lj is an identity matrix, Sj oc (A^)" 1 , so that the uncertainty matrix S^ 1 
is proportional to the background uncertainty matrix A*. For convenience, we label this 
specific form by I-BG, and call c the performance level coefficient (PLC). 

If one ignores Qj in Eq. (125bj) (or if there is no dynamical noise term so that Qj = 0), 
and decomposes A^ and A" as A^ = TiU^Tj and A" = TjU^rf , respectively, where TV, 
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TJ| and U" are suitable matrices, then it can be shown that 



(U 



a\-l 



(33) 



which is equivalent to 



2002 



PhametaL 



1998|) 



;he singular evolutive extended Kalman (SEEK) filter (IHoteit et al 
with a forgetting factor (1 — c). 



b. Case ^LfSiU = c [(A^)" 1 + (Hi)' 1 (Hi)" l/ Hi\ = c (E?)- 1 for < c < 1 



This implies that 



li 



^LjLf )~ 1 L i (S")" 1 (Lf)(LjLf )~ 1 , 



(34) 



and 



a\-l 



a\-l 



(35) 



In particular, if Lj is an identity matrix, S, cx (S") _1 , so that the uncertainty matrix S^" 1 
is proportional to the analysis uncertainty matrix that is updated from through an 
EnKF method. We label this specific form by I- ANA. 



Eq. (1351) means that 



A? 



i™ 



(36) 



which i s equivalent to the coyarian ce inflation technique used in lAnderson and Anderson 



fll999h : 



Whitaker and Hamilll (|2002| ) and other similar works if one lets 1/(1 — c) = (1 +S) 2 , 
with 5 being the inflation factor. Note that, in the EnKF method, one may first update 
the background mean to the analysis counterpart, and then generate an analysis ensemble, 
whose error covariance is (implicitly) inflated. In doing so, covariance inflation does not 
affect the computation of the Kalman gain at the same cycle. Instead, it affects the Kalman 
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gain in the next cycle, since the background covariance at the next cycle will be inflated by 
conducting covariance inflation. In contrast, in the EnTLHF, the Gain matrix is directly 
affected by the analysis uncertainty matrix at the same assimilation cycle. 

c. Case LjSiLi = l mx 
This implies that 

S i = (L i Lf)-\ (37) 

and 

(A?)" 1 = (S?)- 1 - Ti I mi • (38) 

In particular, if Lj is an identity matrix, Sj = I mx , so that the uncertainty matrix S^" 1 is also 
equal to the identity matrix I mx . We label this specific form by I-MTX. 

Suppose that the analysis covariance E" obtained by an EnKF method is decomposed as 

= E?D?(E?) T (39) 

through a singular value decomposition (SVD), where = diag^i, • • ■ , cr^ mx ) is a diagonal 
matrix consisting of the eigenvalues crjj of (j — 1, • • • , m x ), which are arranged in a non- 
ascending order, i.e., a it j > > for j < k, and Ef = [e^, • • • , e ijm J is the matrix 
consisting of the corresponding (orthonormal) eigenvectors e^- (j = 1, • • • ,m x ). Then in 
order to make (A") _1 in Eq. fl3"8"|) positive definite, it is sufficient to let 

<^i-7i>0, (40) 

which means 

7i < — , (41) 
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i.e., 7j is less than the inverse of the maximum eigenvalue of We can also introduce a 
PLC to the filter in this case by writing 7^ = c/cr^x, with < c < 1. 

Given a suitable value of 7^, A" can be evaluated as follows. Using Eq. f|5S|) . one has 

At 1 = (srr 1 - 7^ 

= E?p?)- 1 (E?) T - 7i E?(E?) ;r (42) 
= E»(A^)- 1 (E^) T , 

where Af is a diagonal matrix consisting of the eigenvalues rjij (j — 1, • ■ ■ , mj, with 



J- — liVij 1 — c (Tij/a^x 

Thus one has 

A i » = E i »AJ(E?) T . (44) 

Through Eqs. ( 143]) and ( l44l) . one can see that, the analysis uncertainty matrix A" ob- 
tained by the EnTLHF has a similar structure to that of S" (cf. Eq. (|39|) ). The eigenvalues 
of A" are the inflations of the corresponding eigenvalues a^j of X" if cr^j > 0, or remain 
unchanged if er^- = 0. 



Ottetal 



f l2004h . where the 



Remarks: A similar eigenvalue inflation technique was used in 
authors increased all the eigenvalues of Y^l by a constant e« (relative to the ith cycle), so 
that after inflation, the eigenvalues of A" are equal to <Jij + e», but with the same associated 
eigenvectors as in One may establish the connection of this inflation technique to the 
EnTLHF by setting 7i Lf S;L; = (S^)" 1 - (S^ 1 + ej^)" 1 in Eq. (J3U]). Nominally, in light 
of Eq. fl30|) . one may establish the connection of any inflation technique to the EnTLHF by 
solving the equation 

-,WS,L, = (£< 
24 



,f S.L, = (S^)- 1 - (A")" 1 . 



As one more example, one may consider the scenario "fiLf SjLj = (£") 1 — (A") 1 
ciUifiRi)- 1 ^ for < c < 1. This leads to 



li 



'MM 



- 1 L,(^) T (R,)^ 1 ^(Lf)(L,Ln" 1 



(45) 



and 



a\-l 



b\-l 



[l-c){n i ) T {R i )- 1 'Hi. 



(46) 



To our knowledge, covariance inflation based on Eq. ( 146|) is not used in the literature, possibly 
because it is more natural to conduct covariance inflation through the background or analysis 
ensemble of the system state, rather than through the ensemble of observation. However, 
Eq. (I46p might provid e an alternative po i nt of view to explain the under-performance o. 



the stochastic EnKF (IBurgers et al 



Houtekamer and Mitchell 



1998 



Evensen 



1994 



Evensen and van Leeuwen 



1998) in comparison to the deterministic ones (Anderson 



Bishop et al.ll200ll ; IWhitaker and Hamillll2002l ) in certain situations (e.g. JWhitaker and Hamill 



1996 



2001 



20021 ). In the stochastic EnKF, one generates an ensemble of surrogate observations based on 



the observation distribution. Due to the effect of small sample size, the sample covariance of 
the surrogate observations will under-represent the original covariance of observation. This 
is equivalent to letting 1 — c > 1 in Eq. (j4"6"]) (i.e., c < 0), which implies the choice of 7* < 
in Eq. fl5D|) . As discussed previously, the negativeness of 7, means that, instead of being 
conservative, the filter designer is confident in the estimation accuracy of the background. 
Hence when updating the background to the analysis, more weight is allocated to the back- 
ground, rather than to the observation, which may deteriorate the performance of the filter 



if there exist more uncertainties in the background than in the observation. 
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5. Numerical examples 

We conduct a series of numerical experiments to assess the relative robustness of the 
TLHF/EnTLHF in comparison to the corresponding KF/EnKF method without inflation. 
In all experiments, we estimate the full state vectors so that the transform matrix Lj = I mx . 



a. Experiments with a linear model 

In the first set of experiments, we consider a one-dimensional regression model, governed 
by the equation 

x i+ i = 1 + 0.5x.i - Q.lx? + f(xi] k, h, d) + Ui , (47) 

where the random variable itj follows the Gaussian distribution with mean zero and variance 
1 (denoted by U{ ~ N(ttj : 0, 1)), /(xf, k, h, d) is a discrete carbox function, starting at the 
time index k, with a jump height h and a width d, i.e., 



h, for i — k, k + 1, ■ • • , k + d — 1 



f(xi]k,h,d) = < 



(48) 

0, otherwise , 



In the experiment we let four jumps occur with k = 200, 400, 600 and 800, respectively, 
with the objective to verify that the TLHF can perform (almost) equally well for the jumps 
occurred at different times. For each jump, we fix d = 20, but let the jump height h be 10 
or 30. For illustration, a time series generated by Eq. (147j) with h = 10 is shown in Fig. (TjQ). 
We suppose that the associated observation system is characterized by 

z/i = Xi + Vi, (49) 
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where Vi ~ N(vi : 0, 1). In the experiment, we let the time index i — 1,2, ■■ • , 1000, and 
record the observation every time step. To reduce statistical fluctuations in estimation, we 
conduct the experiment for 100 times, each time with different initial conditions (drawn at 
random), hence different truths and observations. 

For assimilation, we suppose that one uses the following imperfect model 



as the dynamical system, which ignores the nonlinear term — O.lxf and fails to capture the 
abrupt change of regime. Because of the linearity in Eqs. f|49|) and fl50|) . we can apply 
both the KF and the TLHF to assimilate the system. Given the covariance £j (of a scalar 
variable) evaluated by the KF at the ith cycle, in the TLHF one only needs to replace it by 
Aj = Sj/ (1 — 7jSj) = Sj/(1 — c) (with 7^ = c/Sj, c G [0, 1) in the scalar system) and change 
the gain matrix accordingly, while the other steps are the same as those in the KF. Clearly, 
in the scalar system, the specific form I-ANA in §4bl is equivalent to the specific form I-MTX 
in 

We use the average root mean squared error (average RMSE) to measure the perfor- 
mance of the filter. For an m x -dimensional system, the RMSE e; of an estimate Xj = 
(xi t i, ■ ■ ■ ,x ijmx ) T with respect to the state vector x* = (a^i, • • • ,Xi,m x ) T at time instant i is 
defined as 



The average RMSE (RMSE for short) at time instant i over 100 simulations is defined 



We also define the time mean RMSE e as the average of the (average) RMSE over the 



x i+ i = 1 + 0.5xi + Ui 



(50) 




(51) 



as ii = X^'=i e i /100, where e\ denotes the RMSE at time instant % in the jth simulation. 
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time horizon (assimilation time window) [1, ■ • • , N], i.e., e = Yli=i (N = 1000 in this 
experiment). 

Fig.[2]plots the RMSE of the KF over the time horizon [1, 1000] in the case h = 10 (upper 
panel), and that in the case h = 30 (lower panel). In both cases, the KF achieves a relatively 
low RMSE during the period without any abrupt jump. However, when the abrupt jumps 
occur, the RMSE of the KF rises sharply in response. 

Figs. El - E plot the RMSE differences between the TLHF of I-BG with different PLC 
values and the KF. Throughout this work, we use the RMSEs of the KF as the baselines, and 
the RMSE differences are defined as the RMSEs of the TLHF subtracted by the correspond- 
ing ones of the KF. In all these figures, the upper panels correspond to the case h = 10, and 
the lower ones the case h = 30. At c = 0.1 (Fig. |3]), when there is no abrupt jump, the RM- 
SEs of the TLHF and the KF are nearly indistinguishable, so that their RMSE differences 
are almost zero. However, when the abrupt jumps appear, the RMSEs of the TLHF do not 
rise as sharply as those of the KF, so that their RMSE differences become negative, suggest- 
ing that the TLHF has relatively more robust performance than the KF during the abrupt 
jumps. At c = 0.5 (Fig. HJ), the RMSE differences during the periods with the abrupt jumps 
become larger, while those during the periods without any abrupt jump remain close to 
zero. Further increasing c to 0.9 (Fig. EJ), the performance of the TLHF becomes remarkably 
better than the KF during the periods with abrupt jumps. The RMSEs of the TLHF appear 
insensitive to the presence of the abrupt jumps, which is not the case for the KF. However, 
there is also a price for the TLHF to achieve this. During the periods without the abrupt 
jumps, the TLHF performs worse than the KF so that their RMSE differences are slightly 
above zero. Moreover, the divergence of the TLHF is spotted for time indices i > 870. The 
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occurrence of the divergence is possibly due to the fact that the PLC is too large, so that 
l/7i becomes less than the minimum threshold 1/7* defined in Eq. ffT2]) . As discussed in 
£ pbl in such situations there is no guarantee to attain a TLHF solution that satisfies the 
inequality in (fl9|) . Instead, divergence of the filter solution may occur as observed in the 
experiment. 

Figs. El-El show the RMSE differences between the TLHFs of I- ANA and I-MTX with 
three different PLC values (equivalent to each other in scalar systems), and the KF. Similar 
results are observed. At c = 0.1, the RMSEs of the TLHF and the KF are almost indistin- 
guishable when there is no abrupt jump, so that their RMSE differences are very close to 
zero. The TLHF performs better again than the KF when the abrupt jumps occur. At a 
larger PLC value, c = 0.4, the TLHF performs remarkably better than the KF during the 
period of abrupt jumps, but at the cost of slightly worse performance than the KF during 
the period without any abrupt jump. When further increasing c to 0.6, the performance 
of the TLHF deteriorates in comparison with the choice c = 0.4. More investigations (not 
reported here) show that a larger value (c > 0.6) leads to even worse performance. 

To summarize, our experiment results show that, for a relatively small PLC, the KF 
and the TLHF have close performance. This is expected, since the TLHF with c = 
reduces to the KF as we have noted in §31 As c increases, the TLHF exhibits a better 
performance than the KF when there are relatively large uncertainties. However, when 
there only exist relatively small uncertainties in assimilation, a too large c (hence too much 
uncertainty inflation) may also make the TLHF appear over-conservative and deteriorate 
the filter performance (or even diverge). This is because, with relatively small uncertainties, 
the backgrounds also provide useful information, and thus should not be under-weighted. To 
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mitigate this problem, one possible strategy is to use a relatively small value of c to make the 



TLHF less conservative when there only exists relatively small uncertainties, and a larger one 
when there exhibit more uncertainties. This is essentially a strategy that conduct s adaptive 



covariance i n flation, as has already been investigated in som e works, for example 



(12007 



2009|); 



Hoteit et al. 



(|2Q02|); 



Hoteit and Pham 



Anderson 



( 120041 ) . From our earlier discussion in 



§ I5H the adaptive inflation problem can be solved under the framework of the HF with an 
additional optimality criterion (e.g, minimum variance or maximum likelihood), which will 
be investigated in the future. 



b. Experiments with a nonlinear model 



In t 



re second set of experimen ts, we consider the 40-dimensional Lorenz and Emanuel 



model (ILorenz and Emanuel 
given by 

dx 



19981 ) (LE98 model hereafter), whose governing equations are 



— = (x i+1 - Xi- 2 ) Xi-i -Xi+ F(t), i = !,-•• ,40. 



(52) 



The quadratic terms simulate advecti on, the 



F acts as the external forcing term ( ILorenz 



inear term represents internal dissipation, and 



19961 ). For consistency, we define X-\ = 239, 



xq = £40, and £41 = x\. We suppose that the true value of the parameter F is 8 for 
t > 0, but in assimilation one may choose other values for F, which thus yields a potential 
parameter mismatch. In our experiments we consider two scenarios, with F = 6 and F — 8, 
respectively. 

We use the fourth-order Runge-Kutta method to integrate (and discretize) the system 
from time to 250, with a constant integration step of 0.05 (overall 5001 integration steps). 
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We then use the following system 



Yi = Xj + v.; 



(53) 



to record the observation of the state vector Xj = (sc^i, 2^2, ■ ■ • ,£;,4o) T at time instant i, 
where Vj follows the Gaussian distribution iV(vj : 0, 140), with I 40 being the 40-dimensional 
identity matrix. The observations are made for every 4 int egration steps . 



We use the ensemble transform Kalman filter (ETKF) (IBishop et al. 



20011 ) to construct 



the EnTLHF. The ETKFs with I-BG and I-ANA are constructed by inflating the back- 



ground ensembles and the an a 



Anderson and Anderson! (119991 ); 



ysis ensembles, respectively in a similar way to that in 



Whitaker and Hamilll (120021 ) . To construct the ETKF with 



I-MTX, one needs to evaluate the analysis covariances, conduct SVDs, and then inflate the 
associated eigenvalues. In high-dimensional systems, conducting SVDs on the analysis co- 
variances makes the ETKF with I-MTX computationally less efficient than its I-BG and 
I-AN A counterparts. However, it is possib le to implement the I-MTX form in the SEEK 



filter ( IHoteit et all 12002 



Pham et al. 



19981 ) without significant increase of computational 
cost, since in this case all such SVDs can be conducted on the matrices updated by Eq. ()3"3~j) . 
whose dimension is determined by the ensemble size in assimilation. 

In our experiments we let the ensemble size n = 10 and vary the PLC values. To reduce 
statistical fluctuations, for each PLC value c we repeat the experiments for 20 times, each 
time with a randomly drawn initial background ensemble (with 10 members). In practice, 



Hamill et al. 


2009; 


Van Leeuwen 


2009) 



20091 ). Since in our experiments our 



objective is to assess the relative robustness of the EnTLHF, we choose not to conduct 
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covariance localization to avoid complicating the analysis of our experiment results. In what 
follows, we examine the time mean RMSE of the EnTLHF as a function of the PLC value 
c, with c G [0, 0.1, 0, 2, • • • , 0.9]. The ETKF is treated as a special case of the EnTLHF with 
c = 0. 

Fig. [9] plots the time mean RMSEs of the ETKF with I-BG. The result in case of F = 6 
is marked with the dash-dotted line, and that in the case of F = 8 with the dotted one. 
When F — 6, the time mean RMSE appears to be a monotonically decreasing function with 
respect to c. When F — 8, the time mean RMSE tends to decrease until it reaches c = 0.8. 
After that, the time mean RMSE slowly rises. In both cases, all time mean RMSEs with 
c > are lower than that of the ETKF (c = 0). 

Similar results of the ETKF with I- ANA are observed in Fig. [10J For both cases with 
F = 6 and F — 8, their time mean RMSE are monotonically decreasing functions with 
respect to c, and all time mean RMSEs with c > are lower than that of the ETKF (c = 0). 

Fig. [TJJ shows the time mean RMSEs of the ETKF with I-MTX. When F = 6, the time 
mean RMSE decreases monotonically until it reaches c = 0.4. After that, the time mean 
RMSE rises rapidly. Moreover, filter divergence occurs for c > 0.6, possibly for the same 
reason as explained in the previous section. The result of F = 8 is similar: the time mean 
RMSE decreases until c = 0.5, and then increases as c continues growing. Filter divergence 
also occurs when c > 0.6. Compared to the ETKF (c = 0), the time mean RMSEs with 
c > are lower until c reaches the turnaround point. 

Through the above experiments, we have shown that, with suitable PLC values, the 
ETKFs of all three specific forms, namely, I-BG, I-ANA and I-MTX, exhibit relative ro- 
bustness in comparison with the ETKF without any covariance inflation, which is consis- 
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tent with the observations in the literature that an EnKF method with suitable covariance 



inflation may per 



Hamill et al. 



2009 



orm better than th at without any covariance inflation (see, for example, 



Van Leeuwen 



2009 



and the references therein). Different inflation schemes 



may result in different filter performance. For instance, the ETKF with I-ANA appears to 
have better performance than the other two schemes. The validity of this conclusion may 



depend on the system in assimilation, though, and may need to be verified case by case. 



6. Discussion and conclusion 

In this work we considered the applications of the KF and the HF to state estimation 
problem. We discussed the similarity and difference between the KF and the HF, and 
showed that the KF can be considered as a special case of the HF with the performance 
level equal to zero. For convenience of applying the filtering theory to sequential data 
assimilation, we introduced a variant, called the time-local HF, in which we suggested to 
solve the constraints in the HF locally (in time). Analogous to the EnKF methods, we 
proposed the ensemble version of the TLHF, called the ensemble time-local HF (EnTLHF), 
and showed that the EnTLHF can be constructed based on the EnKF. In addition, we 
established the connections of some specific forms of the EnTLHF to some EnKF methods 
equipped with certain covariance inflation techniques. 

Compared to existing works on covariance inflation in the EnKF, the ifoo filtering theory 
provides a theoretical framework that unifies various inflation techniques in the literature, 
and establishes the connection between covariance inflation and robustness. The filtering 
theory also provides an explicit definition of robustness and the associated mathematical 
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description. Conceptually, this leads to the possibility to recast the problem of optimal 
covariance inflation as an optimization problem with multiple objectives, although further 
investigations will be needed for practical considerations. In addition, since the definition of 
robustness is filter- independent, the robustness property may be integrat ed into other types 



of nonlinear filters, e.g., the partic le filter or the Gaussian sum filter ( iHoteit et al. 



Luo et al 



2010 



Van Leeuwen 



2008 



20031 ). by imposing constraints similar to that in Eq. ( Fl9l) . In 



our opinion, it might be less obvious to see how the above extensions can be made from the 
point of view of covariance inflation. 

Through numerical experiments, we verified the relative robustness of three specific forms 
of the TLHF/EnTLHF in comparison with the KF/ETKF without covariance inflation. 
There are also some issues that have not been fully addressed in this work, for instance, 
the optimal choice of the performance level coefficient in conducting uncertainty inflation. 



Further investigations in these aspects will be considered in the future. 
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of Figures 

A time series generated by Eq. ( 1471) . 

RMSE of the KF in assimilating the regression model. The jump heights h 
are 10 (upper panel) and 30 (lower panel), respectively. 

RMSE difference between the TLHF of I-BG and the KF in assimilating the 
regression model. The jump heights h are 10 (upper panel) and 30 (lower 
panel), respectively. In the TLHF, the PLC c = 0.1. 

RMSE difference between the TLHF of I-BG and the KF in assimilating the 
regression model. The jump heights h are 10 (upper panel) and 30 (lower 
panel), respectively. In the TLHF, the PLC c = 0.5. 

RMSE difference between the TLHF of I-BG and the KF in assimilating the 
regression model. The jump heights h are 10 (upper panel) and 30 (lower 
panel), respectively. In the TLHF, the PLC c = 0.9. 

RMSE difference between the TLHF of I-ANA and I-MTX and the KF in 
assimilating the regression model. The jump heights h are 10 (upper panel) 
and 30 (lower panel), respectively. In the TLHF, the PLC c = 0.1. 
RMSE difference between the TLHF of I-ANA and I-MTX and the KF in 
assimilating the regression model. The jump heights h are 10 (upper panel) 
and 30 (lower panel), respectively. In the TLHF, the PLC c = 0.4. 
RMSE difference between the TLHF of I-ANA and I-MTX and the KF in 
assimilating the regression model. The jump heights h are 10 (upper panel) 
and 30 (lower panel), respectively. In the TLHF, the PLC c = 0.6. 
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9 Time mean RMSE of the EnTLHF of I-BG as a function of the PLC in 
assimilating the L96 model. The values of the parameter F are 6 (dash- 
dotted) and 8 (dotted), respectively. The EnTLHF reduces to the ETKF 
when the PLC c = 0. \50 

10 Time mean RMSE of the EnTLHF of I-ANA as a function of the PLC in 
assimilating the L96 model. The values of the parameter F are 6 (dash- 
dotted) and 8 (dotted), respectively. The EnTLHF reduces to the ETKF 
when the PLC c = 0. |si 

11 Time mean RMSE of the EnTLHF of I-MTX as a function of the PLC in 
assimilating the L96 model. The values of the parameter F are 6 (dash- 
dotted) and 8 (dotted), respectively. The EnTLHF reduces to the ETKF 
when the PLC c = 0. Filter divergence occurs for c > 0.6. [52 
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Time step 

Fig. 1. A time series generated by Eq. (1471) . 
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RMSE of the KF with h = 10 
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Fig. 2. RMSE of the KF in assimilating the regression model. The jump heights h are 10 
(upper panel) and 30 (lower panel), respectively. 



43 



RMSE difference between the HF (l-BG) and the KF with h = 10 and c = 0.1 
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Fig. 3. RMSE difference between the TLHF of I-BG and the KF in assimilating the 
regression model. The jump heights h are 10 (upper panel) and 30 (lower panel), respectively. 
In the TLHF, the PLC c = 0.1. 
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RMSE difference between the HF (l-BG) and the KF with h = 10 and c = 0.5 
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Fig. 4. RMSE difference between the TLHF of I-BG and the KF in assimilating the 
regression model. The jump heights h are 10 (upper panel) and 30 (lower panel), respectively. 
In the TLHF, the PLC c = 0.5. 
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RMSE difference between the HF (l-BG) and the KF with h = 10 and c = 0.9 
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Fig. 5. RMSE difference between the TLHF of I-BG and the KF in assimilating the 
regression model. The jump heights h are 10 (upper panel) and 30 (lower panel), respectively. 
In the TLHF, the PLC c = 0.9. 
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RMSE difference between the HF (l-ANAJ-MTX) and the KF with h = 10 and c = 0.1 
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Fig. 6. RMSE difference between the TLHF of I-ANA and I-MTX and the KF in assimi- 
lating the regression model. The jump heights h are 10 (upper panel) and 30 (lower panel), 
respectively. In the TLHF, the PLC c = 0.1. 
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RMSE difference between the HF (l-ANAJ-MTX) and the KF with h = 10 and c = 0.4 
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Fig. 7. RMSE difference between the TLHF of I-ANA and I-MTX and the KF in assimi- 
lating the regression model. The jump heights h are 10 (upper panel) and 30 (lower panel), 
respectively. In the TLHF, the PLC c = 0.4. 
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RMSE difference between the HF (l-ANAJ-MTX) and the KF with h = 10 and c = 0.6 
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Fig. 8. RMSE difference between the TLHF of I-ANA and I-MTX and the KF in assimi- 
lating the regression model. The jump heights h are 10 (upper panel) and 30 (lower panel), 
respectively. In the TLHF, the PLC c = 0.6. 
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Time mean RMSEs of the EnTLHF (l-BG) with different F 
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Fig. 9. Time mean RMSE of the EnTLHF of I-BG as a function of the PLC in assimi- 
lating the L96 model. The values of the parameter F are 6 (dash-dotted) and 8 (dotted), 
respectively. The EnTLHF reduces to the ETKF when the PLC c = 0. 
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Time mean RMSEs of the EnTLHF (l-ANA) with different F 
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Fig. 10. Time mean RMSE of the EnTLHF of I-ANA as a function of the PLC in assim- 
ilating the L96 model. The values of the parameter F are 6 (dash-dotted) and 8 (dotted), 
respectively. The EnTLHF reduces to the ETKF when the PLC c = 0. 
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Time mean RMSEs of the EnTLHF (l-MTX) with different F 
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Fig. 11. Time mean RMSE of the EnTLHF of I-MTX as a function of the PLC in assim- 
ilating the L96 model. The values of the parameter F are 6 (dash-dotted) and 8 (dotted), 
respectively. The EnTLHF reduces to the ETKF when the PLC c = 0. Filter divergence 
occurs for c > 0.6. 
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